// @HEADER
// ***********************************************************************
//
//             SiYuan: A numerical PDE solver
//                 Copyright (2022) YUAN Xi
//
// ***********************************************************************
// @HEADER

#ifndef SCHRODINGERPOTENTIAL_HPP
#define SCHRODINGERPOTENTIAL_HPP

#include "Phalanx_config.hpp"
#include "Phalanx_Evaluator_WithBaseImpl.hpp"
#include "Phalanx_Evaluator_Derived.hpp"
#include "Phalanx_MDField.hpp"

#include "Teuchos_ParameterList.hpp"
#include "Sacado_ParameterAccessor.hpp"
#include "Teuchos_Array.hpp"

#include "Panzer_Dimension.hpp"
#include "Panzer_FieldLibrary.hpp"
#include "Panzer_Evaluator_WithBaseImpl.hpp"

namespace SiYuan {
	
/** 
 * \brief Various Potential Functions used in Shrodinger Equation
 */
template<typename EvalT, typename Traits>
class SchrodingerPotential : public panzer::EvaluatorWithBaseImpl<Traits>,
  	public PHX::EvaluatorDerived<EvalT, Traits>
{
public:
  	typedef typename EvalT::ScalarT ScalarT;

  	SchrodingerPotential(Teuchos::ParameterList& p,
                 const panzer::IntegrationRule&  ir);
  
  	void postRegistrationSetup(typename Traits::SetupData d,
			     PHX::FieldManager<Traits>& vm);
  
  	void evaluateFields(typename Traits::EvalData d);
  
  	//! Function to allow parameters to be exposed for embedded analysis
  	ScalarT& getValue(const std::string &n);

private:

  	ScalarT finiteWallPotential( typename Traits::EvalData workset, const int numDim,
                            const int cell, const int qp );
	ScalarT stringFormulaPotential( typename Traits::EvalData workset, const int numDim,
                            const int cell, const int qp );

	int ir_degree, ir_index;
  	std::size_t numQPs;
  	std::size_t numDims;

  	//! output
    PHX::MDField<ScalarT, panzer::Cell, panzer::IP> m_V; //potential 

  	//! energy parameter of potential, precise meaning dependent on type of potential:
    //  Parabolic case -> confinement energy
    //  Finite Wall case -> barrier height 
  	ScalarT E0;

	//! Origin of coordinate system, used in calculating r
	Teuchos::Array<double> origin;

    //!! specific parameters for string formula
    std::string stringFormula;
    
    //! specific parameters for Finite Wall 
    double barrEffMass; // in [m0]
    double barrWidth;   // in length_unit_in_m
    double wellEffMass;
    double wellWidth; 
    
	  //! constant scaling of potential
  	ScalarT scalingFactor;
  	
  	std::string potentialType;
    std::string potentialStateName;

    //! units
    double energy_unit_in_eV, length_unit_in_m;
};

}

#include "SchrodingerPotential_impl.hpp"

#endif
